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Abstract. We investigate improvements to the algorithm for the computation 
of ideal class groups described by Jacobson in the imaginary quadratic case. 
These improvements rely on the large prime strategy and a new method for 
performing the linear algebra phase. We achieve a significant speed-up and are 
able to compute ideal class groups with discriminants of 110 decimal digits in 
less than a week. 

1. Introduction 

Given a fundamental discriminant A, it is known that the corresponding ideal 
class group C1(A) of the order Oa of discriminant A in K = Q(VA) is a finite 
abelian group that can be decomposed as 

Cl(A)~0Z/diZ, 

i 

where the divisibility condition di\d i+ i holds. In this paper we investigate improve- 
ments in the computation of the group structure of C1(A): that is, determining 
the di, which is of both cryptographic and number theoretic interest. Indeed some 
cryptographic protocols relying on the difficulty of solving the discrete logarithm 
problem (DLP) in imaginary quadratic orders have been proposed [2l[9], and solving 
instances of the DLP is closely related to finding the group structure of C1(A). 

In 1968 Shanks [TO] proposed an algorithm relying on the baby-step giant-step 
method in order to compute the structure of the ideal class group of an imaginary 
quadratic number field in time O (|A| 1 ^ 4+e ), or O (|A| 1 ^ 5+e ) under the extended 
Riemann hypothesis [T3]. This allows us to compute class groups of discriminants 
having up to 20 or 25 decimal digits. Then a subexponential strategy was described 
in 1989 by Hafner and McCurley [SJ. The expected running time of this method is 

e (v / 2+o(l))v / l°g|A|loglog|A|_ 

Buchman and Dullmann [5] computed class groups with discriminants of around 
50 decimal digits using an implementation of this algorithm. An improvement of 
this method was published by Jacobson in 1999 [10 . He achieved a significant 
speed-up by using sieving strategies to generate the matrix of relations. He was 
able to compute the structure of class groups of discriminants having up to 90 
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decimal digits. More recently Sutherland [22] used generic methods in order to 
compute class groups with discriminants having 100 decimal digits. Unlike the 
previous algorithms, this one relies heavily on the particular structure of C1(A) 
thus obtaining variable performances depending on the values of A. 

Our approach is based on that of Jacobson, using new techniques to accelerate 
both the sieving phase and the linear algebra phase; we have obtained the group 
structure of class groups of 110 decimal digit discriminants. 

2. The ideal class group 

In this section we give essential results concerning the ideal class group and 
the subexponential strategies for computing its structure. For a more detailed 
description of the theory of ideal class groups we refer to [5] and [17] ■ In the 
following, A is a non-square integer congruent to or 1 modulo 4, and the quadratic 
order of discriminant A is defined as the Z-modulc 

e> A = z+ A + ^ z. 



We also denote by K the field Q(VA). 

2.1. Description. Elements of C1(A) are obtained from fractional ideals of 0a, 
which arc Z-modules of K of the form: 




where a and b are integers with b = A mod 2 and q is a rational number. The prime 
ideals are the fractional ideals for which there exists a prime number p such that: 

p=pZ + bp + Z or p = pZ (p inert in K). 

Definition 2.1 (Ideal Class group). Let I a be the set of invertible fractional ideals 
of Oa, and Va — {(«) G Ia,® G the subset of principal ideals. We define the 
ideal class group of A as : 

C1(A) :=X A /Va, 

where the group law is the one derived from the multiplication of Z-modules. 

For every a G I a, there exist uniquely determined prime ideals p 1; . . . ,p„ and 
exponents e\, . . . , e„ in Z such that 

a = Pi 1 ■■■Pn- 

Unlike Ta, the ideal class group C1(A) is a finite group. Its order is called the class 
number and usually denoted by h(A). It grows like IAI 1 / 2 " 1 " 6 , as shown in |20j . 

2.2. Computing the group structure. The algorithm for computing the group 
structure of C1(A) is divided into two major phases: relation collection and linear 
algebra. In the first phase, we begin by precomputing a factor base B = {pi, . . . , p„} 
of non-inert prime ideals satisfying A/" (pi) < B, where B is a smoothness bound. 
Then we look for relations of the form 

(«) = pT ■■■p%", 

where a G K. Every n-tuple [ei, . . . , e„] collected becomes a row of what we will 
refer to as the relation matrix A G Z mx ™. We have from [T| the following important 
result: 
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Theorem 2.2. Let A be the lattice spanned by the set of the possible relations. 
Assuming GRH, if B > 6 log 2 A, then we have 

C1(A) ~ Z"/A. 

After the relation collection phase we can test if A has full rank and if its rows 
generate A using methods described in £13.31 If it is not the case then we have to 
compute more relations. From now on we assume that A has full rank and that its 
rows generate A. 

The linear algebra phase consists of computing the Smith Normal Form (SNF) 
of A. Any matrix A in Z nxn with non zero determinant can be written as 



A = V- 1 



( di 


V o 







\ 







£7" 



d n j 



, where for all 1 < % < n and U and V are unimodular matrices in Z nxn . 

The matrix diag(di , . . . , d n ) is called the SNF of A. If m = n and diag(e?i , . . . , d n ) = 
SNF (A) then 

n 

C1(A) ~ 0Z/diZ. 



This reduces the problem of computing the group structure of C1(A) to computing 
the SNF of a relation matrix A in Z nx ™. For an arbitrary A in Z mxn we start by 
computing the Hermite Normal Form (HNF) of A. A matrix H is said to be in 
HNF if it has the shape 



H = 



hi. 







'2,2 



(0) 







\ 



hn,r. 



, where < hij < ha for all j < i and hij — for all j > i. For each matrix A in 
Z mxn there exists a matrix H in HNF and a unimodular matrix W in Z mxm such 
that 

H = WA. 

The upper block of H is a n x n relation matrix whose SNF provides us the group 
structure of C1(A). There is an index I such that h^i = 1 for every i > I. The upper 
left I x / submatrix of H is called the essential part of H . In order to compute the 
group structure of C1(A) it suffices to compute the SNF of the essential part of H, 
which happens to have small dimension in our context. 
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2.3. The use of sieving for computing the relation matrix. The use of 
sieving to create the relation matrix was first described by Jacobson |10j . Here we 
follow the approach of [TT] Chap. 13, which relies on the following lemma: 

Lemma 2.3. If a = (aZ + ^Az) with a > 0, then for all x,y in Z there exists 
b G Ia such that ab € V and 

b 2 - A 

Af(b) = ax 2 + bxy H y 2 . 

The strategy for finding relations is the following: We start with 




whose norm is i3-smooth. Then we choose a sieve radius R satisfying R « v/|A|/2/W(a) 
and we look for values of x £ [-R, R] such that (p(x, 1) is i?-smooth where 

/ \ 2 . b 2 -A 2 

(p(x,y) = ax +bxy-\ y , 

4a 

which allows us to find b = FJ ■ satisfying ab = (7) for some 7 in K. The pi and 
fi are deduced from the decomposition tp(x, 1) = Oi^i*- For more details we refer 
to [IT], Chap 13. This method yields the relation 

(7)=n»r +/i - 

i 

Now given a binary quadratic form <p(x, y) = ax 2 + bxy + cy 2 of discriminant A, 
we are interested in finding values of x £ [-R, R] such that ip(x, 1) is _B-smooth. 
This can be done trivially by testing all the possible values of x, but there is a 
well-known method for pre-selecting some values of x in [-R, R] that are going to 
be tested, namely the quadratic sieve (introduced by Pomerance [E]). It consists 
in initializing to an array S of length 2R + 1 and precomputing the roots r\ and 

r", or the double root of <p(x, 1) mod pi for each Pi < B such that (^j^j 7^ — 1 • 
Then for each x in [— R, R] of the form x — i\ + kpi for some fc, we add |log£>iJ to 
S[x\. At the end of this procedure, if <p(x, 1) is B-smooth, then S[x] ~ \og(p(x, 1). 
As ip(x, 1) ~ y/A/2R, we set a bound 

(1) F = logl J^R I -TlogCp,, 




where T is a number representing the tolerance to rounding errors due to integer 
approximations. We then perform a trial division test on every ip(x, 1) such that 
Six] > F. 



3. Practical improvements 

In this section we describe the improvements that allowed us to achieve a signif- 
icant speed-up with respect to the existing algorithm and the computation of class 
group structures of large discriminants. Our contribution is to take advantage of 
the large prime variants, of an algorithm due to Vollmer [23] for the SNF which had 
not been implemented in the past, and of special Gaussian elimination techniques. 
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3.1. Large prime variants. The large prime variants were developed in the con- 
text of integer factorization to speed up the relation collection phase in both the 
quadratic sieve and the number field sieve. Jacobson considered analogous variants 
for class group computation [10] , but the speed-up of the relation collection phase 
was achieved at the price of such a slow-down of the linear algebra that it did not 
significantly improve the overall time. The main idea is the following: We define the 
"small primes" to be the prime ideals in the factor base and the small prime bound 
as the corresponding bound B± = B. Then we define a large prime bound B 2 - 
During the relation collection phase we choose not to restrict ourselves to relations 
only involving primes p in B but we also keep relations of the form 

(a) =p!...p n p and (a) = pi . . . p n pp' 

for pi in B, and for p,p' of norm less than B 2 - We will respectively refer to them as 

1- partial relations and 2-partial relations. Keeping partial relations only involving 
one large prime is the single large prime variant, whereas keeping two of them is 
the double large prime variant which was first described by Lenstra and Manasse 
[12] . In this paper we do not consider the case of more large primes, but it is a 
possibility that has been studied in the context of factorization [T3] . 

Partial relations may be identified as follows. Let m be the residue of <p(x, 1) 
after the division by all primes p < B\ , and assume that B 2 < B\ . If m = 1 then 
we have a full relation. If m < B 2 then we have a 1-partial relation. We can see 
here that detecting 1-partial relations is almost for free. If we also intend to collect 

2- partial relations then we have to consider the following possibilities: 

1. m > B\\ 

2. m is prime and m > B%\ 

3. m < B 2 ; 

4. m is composite and B\ < m < B 2 . 

In Cases 1 and 2 we discard the relation. In Case 3 we have a 1-partial relation, and 
in Case 4 we have m = pp' where p = Af(p) and p' = Af(p'). After testing if we are 
in Cases 1, 2, or 3 we have to factorize the residue. We have done that using Milan's 
implementation of the SQUFOF algorithm [TS] based on the theoretical work of [7]. 

Even though we might have to factor the residue, collecting a partial relation is 
much faster than collecting a full relation because the probability that Af(b) is B 2 - 
smooth is much greater than the probability that it is Bi-smooth. This improvement 
in the speed of the relation collection phase comes at a price: The number of columns 
in the relation matrix is much greater, thus preventing us from running the linear 
algebra phase directly on the resulting relation matrix and forcing us to find many 
more relations since we have to produce a full rank matrix. We will see in !j3.2l 
how to reduce the dimensions of the relation matrix using Gaussian elimination 
techniques and in Sj4] how to optimize the parameters to make the creation of the 
relation matrix faster, even though there are many more relations to be found. 

3.2. GAUSSIAN elimination techniques. Traditionally rows were recombined to 
give full relations as follows: In the case of 1-partial relations, any pair of relations 
involving the same large prime p were recombined into a full relation. In the case 
of 2-partial relations, Lenstra [12] described the construction of a graph whose 
vertices were the relations and whose edges linked vertices having one large prime in 
common. Finding independent cycles in this graph allows us to find recombinations 
of partial relations into full relations. 
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In this paper we rather follow the approach of Cavallar [4], developed for the 
number field sieve, which uses Gaussian elimination on columns without distin- 
guishing those corresponding to the large primes from the others. One of the main 
differences between our relation matrices and the matrices produced in the number 
field sieve is that our entries are in Z rather than F2, thus obliging us to monitor the 
evolution of the size of the coefficients. Indeed, eliminating columns at the price of 
an explosion of the size of the coefficients can be counter-productive in preparation 
for the HNF algorithm. 

In what follows we will use a few standard definitions that we briefly recall here. 
First, subtracting two rows is called merging. This is because rows are stored as 
lists of the non-zero entries sorted with respect to the corresponding columns and 
subtracting them corresponds to merging the two sorted lists. If two rows r\ and 
ri share the same prime p with coefficients C\ and C2 respectively then multipling 
T\ by C2 and ri by c\ and merging is called 'pivoting. Finally, finding a sequence 
of pivots leading to the elimination of a column of Hamming weight k is a fc-way 
merge. 

We aim to reduce the dimension of the relation matrix by performing fc-way 
merges on the columns of weight k = 1, . . . , w in increasing order for a certain 
bound w. Unfortunately, the density of the rows and the size of the coefficients in- 
crease during the course of the algorithm, thus obliging us to use optimized pivoting 
strategies. In what follows we describe an algorithm performing A:- way merges to 
minimize the growth of both the density and the size of the coefficients. 

First we have to define a cost function defined over the set of the rows encapsu- 
lating the difficulty induced for the HNF algorithm. In factorization, we want to 
find a vector in the kernel of the relation matrix which is defined over F2 ; the only 
property of the row that really matters is its Hamming weight. In our context, we 
need to minimize the Hamming weight of the row, but we also have to take into 
account the size of the coefficients. Different cost functions lead to different elimina- 
tion strategies. Our cost function was determined empirically: We took the number 
of non-zero entries, counting c times those whose absolute value was above a bound 
Q, where c is a positive number. If r = [ei, . . . , e n ] corresponds to (a) = YiiPT 
then 

C(r)= Yl 1 + c E L 

l<|e<|<Q |ej|>Q 

Indeed as we will see, matrices with small entries are better suited for the HNF 
algorithm described in £13.31 Let us assume now that we are to perform a fc-way 
merge on a given column. We construct a complete graph Q of size fc as follows: 

• The vertices are the rows 

• Every edge linking n and Tj is labeled by C(rfj), where is obtained by 
pivoting Ti and rj . 

Finding the best sequence of pivots with respect to the cost function c we chose 
is equivalent to finding the minimum spanning tree T of Q, and then recombining 
every row r with its parent starting with the leaves of T. 

Unfortunately, some coefficients might grow during the course of column elimi- 
nations despite the use of this strategy. Once a big coefficient is created in a given 
row r, it is likely to spread to other rows once r is involved in another column 
elimination. We must therefore discard such rows as quickly as possible. In our 
implementation we chose to do it regularly: Once we have performed all the fc-way 
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merges for k < 10 • i and i = 1, . . . , w/10 we discard a fixed number K of the rows 
containing the largest coefficients. 

We show in TableQ]the effect of the use of a cost function taking into account the 
size of the coefficients and the regular discard of the worst rows for A = — 4(10 70 + f ) 
with c = 100, Q = 8 and K — 10. We kept track of the evolution of the dimensions 
of the matrix, the average Hamming weight of the rows, and the maximum and 
minimum size of the coefficients. In the first case we use the traditional cost function 
that only takes into account the Hamming weight of the rows and we keep deleting 
the worst rows regularly; this corresponds to taking c = 1 and K = 10. In the 
second case, we use the cost function described above but without row elimination 
by setting c = 100 and K = 0. In the third case, we combine the two (c = 100 and 
K = 10). We clearly see that the coefficients are properly monitored only in the 
latter case. Indeed using a cost function that does not take into account the size of 
the coefficients and just discarding the worst rows regularly seems more efficient in 
terms of reduction of the matrix dimension, but the row corresponding to i = 12 
(that is to say after all the 120-way merges) clearly shows that we run the risk of 
an explosion of the coefficients. 



Figure 1 . Comparative table of elimination strategies 



Without score depending on the size of the coefficients 


i 


Row Nb 


Col Nb 


Average weight 


max 


min 





38752 


45975 


22 


10 


-10 


2 


2334 


1668 


76 


21 


-20 


4 


2123 


1477 


117 


52 


-56 


6 


2028 


1402 


146 


59 


-62 


8 


1951 


1345 


175 


72 


-65 


10 


1890 


1304 


203 


193 


-196 


12 


1836 


1270 


219 


212 


-2147483648 


Without row elimination 


i 


Row Nb 


Col Nb 


Average weight 


max 


min 





38752 


45975 


22 


10 


-10 


2 


2373 


1687 


79 


30 


-40 


4 


2224 


1538 


118 


67 


-50 


6 


2158 


1472 


148 


71 


-132 


8 


2117 


1431 


179 


2648 


-10568 


10 


2097 


1411 


196 


347136 


-337920 


12 


2080 


1394 


214 


268763136 


-173162496 


With adapted score and row elimination 


i 


Row Nb 


Col Nb 


Average weight 


max 


min 





38752 


45975 


22 


10 


-10 


2 


2357 


1691 


76 


17 


-17 


4 


2176 


1530 


114 


27 


-30 


6 


2074 


1448 


149 


37 


-37 


8 


2013 


1407 


177 


43 


-43 


10 


1958 


1372 


199 


44 


-45 


12 


1908 


1342 


224 


54 


-53 



3.3. Vollmer's algorithm for computing the HNF. In [TO] it has been ob- 
served that the algorithm used to compute the HNF of the relation matrix relied 
heavily on the sparsity of the matrix. While recombinations of the kind described 
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in 12 or the techniques of §3.21 reduce the dimensions of the matrix, they also dra- 
matically increase the density of the matrix, thus slowing down the computation 
of the HNF. We had to find an HNF algorithm whose features were adapted to 
our situation. Vollmer described in |23j an algorithm of polynomial complexity de- 
pending on the capacity to solve diophantine linear systems, but not on the density 
of the matrix. It was not implemented at the time because there was no efficient 
diophantine linear system solver available. We implemented Vollmer's algorithm 
using the IML [21 1 library provided by Storjohann. 

Here we give a brief description of the algorithm (for more details we refer to 
[23]). We assume we have created an m x n relation matrix A of full rank. For each 
i < n, we define two matrices 



-4, 



O-mA 



\ 



Ol . 



I \ 



and ei 







V "" J V 1 / 

For each i, we define hi to be the minimal denominator of a rational solution of the 
system 

AiX — 6^, 

this is computed using the function MinCertif iedSol of IML, which is an imple- 
mentation of (Special)MinimalSolution from [16) . and used in [23] for the complexity 
analysis. In [23j it is shown that 



h(A) = Y[hi 



Fortunately, analytic formulae allow us to compute a bound /i* such that 

h* < h(A) < 2h*, 

so we do not have to compute hi for every i £ [1,^]- In addition, the matrices 
produced for the computation of the group structure of C1(A) have small essential 
part, which keeps the number of diophantine systems to solve small (about the same 
size as the number of columns of the essential part) as shown in |23j . 



Algorithm 1 Computation of the class number 

Input: A, relation matrix A of full rank and h* 
Output: h(A) 

h<-l 

i <— n 

while h < hit do 

Compute the minimal denominator hi of a solution of A4 ■ x = 

h <— h- hi 

i <- i — 1 
end while 
return h 



We can compute the essential part of the HNF of A with a little extra effort 
involving only modular reductions of coefficients; we refer to |23j for more details. 
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This part of the algorithm is highly dependent on the performance of the diophantine 
solver we use, which in turn is mostly influenced by the number of columns of the 
matrix and the size of the coefficients. The bechmarks available [3T] show that the 
algorithm runs much faster on matrices with 3-bit coefficients, which is why we took 
coefficient size into account in the cost function for the Gaussian elimination. 

4. Optimization of the parameters 

In this section we proceed to optimize the parameters involved in the relation 
collection phase. Each parameter has an effect on the overall time taken to compute 
the group structure of C1(A). Recall ([TJ giving the bound F; when we collect partial 
relations it should be adapted in the following way: 

F = log (jfjRj -TlogB 2 , 
where B2 is the large prime bound. 

4.1. Optimization of T. The parameter T represents the tolerance to rounding 
errors in the traditional sieving algorithms. Its value is empirically determined, and 
usually lies in the interval [1,2]. In the large prime variant it also encapsulates the 
number of large primes we want to allow. Indeed, if there were no rounding errors 
one would expect this value to be 1 for one large prime and 2 for two large primes. 
In practice, we can exhibit an optimum value which differs slightly from what we 
would expect. In figure [5] we show the overall running time of the algorithm when 
the parameter T varies between 1.5 and 3.5 for the discriminant A = — 4(10 75 + 1). 
The size of the factor base taken is 3250, the ratio BijB\ equals 120, and we allow 
two large primes. 



Figure 2 . Optimum value of T 



.1 15000 



5000 



2.5 
T 



One of the main issues for determining the optimal value of T is that it tends to 
shift when one modifies the value of B\, the rest being unchanged. Indeed, if for 
example B2/B1 = 120 then 

F = log (J^r) -T log 1205!, 
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so when we increase B\ we have to lower T to compensate. Figure [3] illustrates this 
phenomenon on the example A = — 4(10 75 + 1), with two large primes. 



Figure 3. Effect of \B\ on the optimal value of T 



3.5 





+ 
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+ 












+ 

+ 



1.5 

2500 3000 3500 4000 4500 

FB 

In Figure 2] we study the evolution of the optimal value of T for the single and 
double large prime variants on discriminants of the form —4(10" + 1) where n ranges 
between 60 and 80. It appears that, as we expected, the optimal value for the double 
large prime variant is greater than the one corresponding to the single large prime 
variant. This value is between 2 and 2.3 for one large prime and around 2.7 when 
we allow two large primes. 



Figure 4. Optimal value of T when n varies 



2.8 
2.6 - 
2.4 - 
2.2 



□ □ □ □ 



+ + 



2 - + + + 1 large prime + 

2 large primes □ 



1.8 



55 60 65 70 75 80 85 



4.2. The size of the factor base. The optimal size of the factor base reflects 
the trade-off between the time spent on the relation collection phase and on the 
linear algebra phase. This optimum is usually not the size that minimizes the 
time spent on the relation collection phase. To illustrate this, Figure [5] shows the 
time taken by the algorithm for A = -4(10 75 + 1) with B 2 /Bx = 120 and the 
corresponding optimal T. 
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Figure 5. Optimal value of \B\ 
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The optimal size of the factor base increases with the size of the discriminant. 
Figure |6] shows the optimal size of the factor base for discriminants of the form 
—4(10™ + 1) as n ranges between 60 and 80 for both one large prime and two large 
primes. We notice that the single large prime variant requires smaller factor bases 
than without large primes, and bigger factor bases than the double large prime 
variant. 



Figure 6. Optimal value of \B\ when n varies 
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4.3. The ratio B 2 /Bi. Theoretically B 2 should not exceed B\ . In practice, when 
the ratio B 2 jB\ is too high we lose time taking into account partial relations in- 
volving primes that are so large that they are very unlikely to occur twice and to 
lead to a recombination. This phenomenon is known in the context of factorization, 
and 120 is a common choice of value of B2/B1 (see [5]). We ran experiments using 
12, 120 and 1200 as values for the ratio B2/B1. Figure [7] shows the results for 
A = — 4(10 75 + 1) with two large primes. We give the optimum timings for each 
value of the size of the factor base, and compare those values for the three different 
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ratios. It appears that 120 is indeed the best choice, but the performance of the 
algorithm is not highly dependent on this parameter. 



Figure 7. Comparative table of B 2 /Bi 
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5. Computational results 

5.1. Comparative timings. In Figure [5] we give comparative timings in seconds 
between no large primes and the large prime variants for discriminants of the form 
-4(10™ + 1), for n between 60 and 80. We used 2.4GHz Opterons with 16GB of 
memory, and the NTL library with GMP. It appears that we achieved a significant 
speed-up by using the large prime strategy. Direct comparison with previous meth- 
ods based on sieving is hard since the timings available in |10j were obtained on 
296 MHz UtraSPARC-II processors; therefore we just quote that the computation 
of the group structure corresponding to A = — 4(10 80 + 1) took 5.37 days (463968 
CPU time) at the time. We also notice that the double large prime variant does not 
provide an impressive improvement on the overall time for the sizes of discriminant 
involved. The performance is comparable for discriminants of 60 decimal digits and 
starts showing an improvement when we reach 75 digit discrimants. 

Figure 8. Comparative table of the performances (CPU time) 



n 


Large primes 


1 Large prime 


2 Large primes 


60 


374 


284 


280 


65 


1019 


756 


776 


68 


2010 


1489 


1122 


70 


2148 


1663 


1680 


75 


8409 


6669 


5347 


80 


21215 


17123 


14664 



5.2. Large discriminants. In the imaginary case, the largest class groups that 
had been computed using relation collection methods had 90 digits; some 100 dec- 
imal digit discriminant class group structures could be computed using the tech- 
niques of |22) . With the techniques described in this paper, we achieved the com- 
putation of a class group with a 110 decimal digit discriminant. We used 100 Core2 
Duo 2.4GHz Pentium IV processors with 2 GB of memory each for the sieving, 
and one 2.66 GHz Opteron with 64 GB of memory for the linear algebra, which 
is the real bottleneck of this algorithm. Indeed, the sieving phase can be trivially 
parallelized for as many processors as we have and does not require much memory, 
whereas the linear algebra can only be parallelized into the number of factors of 
h that we get from Vollmer's algorithm (around 10 in our examples) and requires 
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a lot of memory. Indeed the limit in terms of matrix dimensions for the diophan- 
tine solver on a 64GB memory computer seems to be around 10000 columns. For 
comparison, in the case of the 110 decimal digit discriminant we had to handle an 
8000-column matrix (after the Gaussian reduction). 



Figure 9. Decomposition of C1(A) for A = -4(10™ + 1) 



n 


decomposition 


100 
110 


C(2) Y x C(1462491779472195274571694315857495335176880879072) 
C(2) n x (7(8576403641950292891121955131452148838284294200071440) 
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